#include "cem/assembler.h"
#include "cem/electromagnetics_system.h"
#include "cem/material_database.h"
#include "cem/mesh_manager.h"
#include "cem/parameter_reader.h"
#include "cem/boundary_condition.h"

#include <gtest/gtest.h>

#include <libmesh/libmesh.h>
#include <libmesh/numeric_vector.h>
#include <libmesh/sparse_matrix.h>
#include <libmesh/linear_implicit_system.h>

#include <random>
#include <string>
#include <fstream>

using cem::Assembler;
using cem::BoundaryConditionManager;
using cem::ElectromagneticsSystem;
using cem::MaterialDatabase;
using cem::MeshManager;
using cem::ParameterReader;

static libMesh::LibMeshInit * g_libmesh_init = nullptr;

// 函数: build_mesh_path
// 作用: 构造 example/case1/3D.msh 的绝对路径.
// 参数: 无
// 返回: 绝对路径字符串
static std::string build_mesh_path()
{
#ifdef PROJECT_SOURCE_DIR
    std::string base = PROJECT_SOURCE_DIR;
    if (!base.empty() && base.back() == '/')
    {
        base.pop_back();
    }
    return base + "/example/case1/3D.msh";
#else
    return std::string("./example/case1/3D.msh");
#endif
}

// 函数: build_temp_json_path
// 作用: 构造当前工作目录下的临时 JSON 文件路径.
// 参数: fname - 文件名
// 返回: 拼接后的路径
static std::string build_temp_json_path(const std::string & fname)
{
    return std::string("./") + fname;
}

// 函数: write_text_file
// 作用: 将给定文本内容写入到指定路径, 覆盖已存在内容.
// 参数: path - 文件路径; text - 文件内容
// 返回: 无 (断言保证写入成功)
static void write_text_file(const std::string & path, const std::string & text)
{
    std::ofstream ofs(path, std::ios::binary);
    ASSERT_TRUE(ofs.is_open()) << "failed to open file for writing: " << path;
    ofs << text;
    ofs.close();
}

// 函数: make_config_text
// 作用: 生成仅欧拉角不同(加 360°)的两份各向异性配置 JSON 文本.
// 参数: add_full_turn - 若为 true, 则在每个角度上加 360 度; 否则使用基准角度.
// 返回: JSON 文本字符串
static std::string make_config_text(bool add_full_turn)
{
    const double s = 30.0 + (add_full_turn ? 360.0 : 0.0);
    const double d = 40.0 + (add_full_turn ? 360.0 : 0.0);
    const double l = 50.0 + (add_full_turn ? 360.0 : 0.0);

    // 说明:
    // - FREQUENCY: 提供一个频率以初始化 omega
    // - PHYSICAL_AREA: 定义两个区域(与示例网格一致), 各向同性基础参数
    // - Anisotropicc: 针对区域 2 指定电导率各向异性, 角度采用 Rz*Rx*Ry 外转序列
    //   由于三角函数周期性, (θ+360°) 应与 θ 等效, 因此装配矩阵应一致.
    std::string json = R"JSON({
        "FREQUENCY": [{"value": 10.0}],
        "PHYSICAL_AREA": {
            "EPSL": [
                {"area": 1, "value": 1.00059},
                {"area": 2, "value": 4.0}
            ],
            "COND": [
                {"area": 1, "value": 1e-10},
                {"area": 2, "value": 1e-2}
            ],
            "MU": [
                {"area": 1, "value": 1.0},
                {"area": 2, "value": 1.0}
            ]
        },
        "Anisotropicc": [
            {"area": 2, "sigma1": 0.02, "sigma2": 0.01, "sigma3": 0.005,
             "strike": %STRIKE%, "dip": %DIP%, "slant": %SLANT%}
        ]
    })JSON";

    auto replace_all = [](std::string & s, const std::string & from, const std::string & to) {
        size_t pos = 0;
        while ((pos = s.find(from, pos)) != std::string::npos)
        {
            s.replace(pos, from.length(), to);
            pos += to.length();
        }
    };

    replace_all(json, "%STRIKE%", std::to_string(s));
    replace_all(json, "%DIP%", std::to_string(d));
    replace_all(json, "%SLANT%", std::to_string(l));
    return json;
}

// 函数: assemble_system_from_json
// 作用: 基于给定 JSON 文本构造系统、载入网格并完成装配.
// 参数: json_text - JSON 配置文本; out_system - 输出系统对象以便后续比较; out_mesh_mgr - 输出网格管理器以保持网格生命周期
// 返回: 无(断言失败将导致用例失败)
static void assemble_system_from_json(const std::string & json_text,
                                      ElectromagneticsSystem & out_system,
                                      MeshManager & out_mesh_mgr)
{
    ASSERT_NE(g_libmesh_init, nullptr);

    // 载入网格
    MeshManager mesh(g_libmesh_init->comm());
    mesh.load_gmsh_mesh(build_mesh_path());

    // 将 JSON 文本落盘为临时文件并加载
    const std::string tmp_path = build_temp_json_path("tmp_rotation_config.json");
    write_text_file(tmp_path, json_text);

    ParameterReader reader;
    ASSERT_TRUE(reader.load_from_file(tmp_path));

    // 初始化系统
    MaterialDatabase material_db;
    BoundaryConditionManager bc_mgr; // 空边界(装配内会自动添加锚定以避免奇异)
    out_system.initialize(reader, mesh, material_db, bc_mgr);

    // 装配
    Assembler assembler;
    assembler.assemble(out_system);

    // 维持网格生命周期
    out_mesh_mgr = std::move(mesh);
}

// 用例: 旋转一致性 — 欧拉角增加 360° 不改变装配矩阵在线性算子作用下的结果
TEST(RotationConsistencyTest, matrix_action_invariance_under_full_turn)
{
    ElectromagneticsSystem system_base;
    ElectromagneticsSystem system_turn;
    MeshManager mesh_base(g_libmesh_init->comm());
    MeshManager mesh_turn(g_libmesh_init->comm());

    assemble_system_from_json(make_config_text(false), system_base, mesh_base);
    assemble_system_from_json(make_config_text(true), system_turn, mesh_turn);

    auto & sys1 = system_base.linear_system();
    auto & sys2 = system_turn.linear_system();

    // 矩阵关闭以确保可用于乘法
    sys1.matrix->close();
    sys2.matrix->close();

    // 构造若干确定性随机向量, 检查 A1*x 与 A2*x 的相对误差
    const libMesh::Parallel::Communicator & comm = g_libmesh_init->comm();

    using number_t = libMesh::Number;
    std::unique_ptr<libMesh::NumericVector<number_t>> x;
    std::unique_ptr<libMesh::NumericVector<number_t>> y1;
    std::unique_ptr<libMesh::NumericVector<number_t>> y2;

    x = libMesh::NumericVector<number_t>::build(comm);
    y1 = libMesh::NumericVector<number_t>::build(comm);
    y2 = libMesh::NumericVector<number_t>::build(comm);

    const libMesh::dof_id_type n = sys1.rhs->size();
    ASSERT_EQ(n, sys2.rhs->size());

    x->init(n, false, libMesh::PARALLEL);
    y1->init(n, false, libMesh::PARALLEL);
    y2->init(n, false, libMesh::PARALLEL);

    std::mt19937 rng(42);
    std::uniform_real_distribution<double> dist(-1.0, 1.0);

    const int trials = 5;
    for (int t = 0; t < trials; ++t)
    {
        // 填充 x（仅写入本地拥有区间）
        x->zero();
        const auto first = x->first_local_index();
        const auto last  = x->last_local_index();
        for (libMesh::dof_id_type i = first; i < last; ++i)
        {
            number_t v = number_t(dist(rng));
            x->set(i, v);
        }
        x->close();

        // y1 = A1 * x, y2 = A2 * x
        sys1.matrix->vector_mult(*y1, *x);
        sys2.matrix->vector_mult(*y2, *x);

        // 计算差的相对范数
        y1->add(static_cast<number_t>(-1.0), *y2); // y1 <- y1 - y2
        const double diff_norm = y1->l2_norm();

        // 恢复 y1 原值以避免影响下一次试验(可选)
        y1->add(static_cast<number_t>(1.0), *y2);

        const double n1 = y1->l2_norm();
        const double n2 = y2->l2_norm();
        const double denom = std::max({1e-14, n1, n2});
        const double rel = diff_norm / denom;

        EXPECT_LT(rel, 1e-10) << "trial " << t << ": relative difference = " << rel;
    }
}

int main(int argc, char ** argv)
{
    ::testing::InitGoogleTest(&argc, argv);
    libMesh::LibMeshInit init(argc, argv);
    g_libmesh_init = &init;
    return RUN_ALL_TESTS();
}